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ABSTRACT 

A  stable   numerical   procedure   for   solving   the   coupled, 
nonlinear  equations   of   electron   and   ion   conservation    and 
electrical    potential    (Poisson's   equation),    is   described.      A 
two-dimensional   model   with    periodic   active    sites   on    a   flat 
plate   is  utilized   to   obtain    both  qualitative  and   quantita- 
tive  results  which  clearly    illustrate   the   self-generating 
sheath    and   ambipolar    regions  adjacent    to   a   non-emitting 
electrode.      The   active  sites  are  characterized   as    voltage 
sources    and   by  electron  densities   depressed   from   both   the 
non-active   wall   and  the  free-stream   values.      Various   cases 
of    species   density  at   the   nonactive   wall   are  evaluated. 
Recombination/ion izat ion   plays    an   important   role   in    estab- 
lishing  the  boundary    layer   behavior   in  the   ambipolar    region 
and   in    the    dimensionality   of  the   problem    formulation. 
Application   is   for   Nitrogen    gas  at   one   amagat  in    a   pulsed 
discharge  of  about   50   micro-seconds  in  duration. 

Results  clearly   demonstrate  the   boundary  layer   nature  of 
the   species   density   at  the   electrode.      A    sheath   length  of   38 
microns   about  a    35-volt  active   site  and  5-20   microns    along 
the   non-active   portion  of   the  electrode   is   established. 
Joule   heating  is   determined    to   be   not   important. 
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I.    INTRODUCTION 

A.       BACKGROUND 

Losses    at   the  electrode    boundary   regions  play   an    impor- 
tant  role   in   high   energy   density   molecular   gas   flow   lasers. 
Of   particular  interest  to   this   study   is  the   Carbon-Dioxide 
glow-discharge   plasma,    of  which   Nitrogen   is   the   principal 
energy    bearing   constituent. 

The   main  energy    loss   mechanisms   in  a   collision   dominated 
plasma    are    found   in   the  sheath    and  ambipolar   regions.      The 
sheath   voltage   loss    occurs    directly   as  a    result   of   Debye 
shielding,    which    creates  a    space  charge   layer  abutting   an 
electrode  or   wall.      The   ambipolar   losses    are  associated    with 
ambipolar   diffusion,    relating   to  the    charged  species    density 
being   less   than   the   free-stream   or   equilibrium   density. 
This   loss   exists    even   in  the  absence   of   current   flow. 

An    initial   premise  of   probe   theories    that  takes  account 
of   collisions  between   charged  particles,    atoms   and   molec- 
ules,   is   that  the    plasma   is    guasi-neutral    up   to   a    certain 
distance   from   the    probe  and   that   the   electric  fields    are 
localized   within    a   narrow    sheath  region  adjacent   to   the   sur- 
face  of    the   probe.      The  problem   of   making   the  transition 
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from  the  quasi-neutral  solution  to  the  sheath  solution  was 
usually  resolved  through  matching  techniques.   The  descrip- 
tion of  a  collision  dominated  sheath  together  with  the  qua- 
si-neutral ambipolar  region  presents  an  imposing  set  of  cou- 
pled nonlinear  partial  differential  equations,  for  which 
there  exists  no  known  exact  solution.   This  study  reports  a 
numerical  method  that  describes  the  sheath  and  ambipolar 
regions  in  a  self-consistent  fashion. 

In  the  hydrodynamic  approximation  (i.e.,  collision  domi- 
nated), the  dimensions  of  the  space  charge  region  "/^j"  (the 
sheath  length)  is  larger  than  the  mean  free  path  "  /  "  of  the 
species  particles.   Additionally  the  sheath  characteristic 
length  is  very  much  smaller  than  the  characteristic  length 
of  the  plasma  boundary  layer.   This  puts  our  study  in  the 
collisionally  thin  sheath  region,  [Ref.  1  ]. 

However,  the  plasma  is  not  considered  to  be  in  thermal 
equilibrium.   In  fact  for  plasmas  in  thermal  nonequilibrium, 
the  electron  temperature  is  governed  by  the  electron  energy 
equation,  which  includes  the  effects  of  the  electron  energy- 
conduction  and  electric  field  as  well  as  collisional  energy 
exchange  between  electrons  and  heavy-gas  particles.   The 
consideration  of  negligible  energy  exchange  during  colli- 
sions (thermally  frozen)  does  not  generally  result  in  a 
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constant  electron  temperature  (T€) ;  particularly  in  the 
sheath  region,  even  if  other  properties  such  as  the  ion 
temperature  (T( )  are  constant.   Therefore  the  electron 
energy  equation  is  necessary  even  in  the  case  of  the  ther- 
mally frozen  situation. 

An  ionizer/sustainer  discharge  pumped  laser  requires 
that  relatively  high  E/n  values  be  distributed  as  uniformly 
as  possible  throughout  the  bulk  of  the  gas  medium.   The 
electrode  regions  however  represent  a  substantially  higher 
and  more  inhomogeneous  E/n  which  may  impose  an  instability 
on  the  remainder  of  the  discharge.   In  the  region  of  the 
electrode,  the  E/n  is  significantly  higher  than  the  undis- 
turbed plasma.   In  this  region  the  elevated  electron  temper- 
ature, which  is  a  direct  function  of  the  sustainer  field, 
will  increase  ionization,  which  could  lead  to  arcing  and 
breakdown. 

In  the  so-called  vacuum  arc  [Ref.  2],  anode  spots  appear 
under  certain  operating  conditions.   These  anode  spots  are 
of  practical  importance,  having  been  also  observed  in  high 
density  plasmas  "Ref.  3].   They  are  typically  micron  or 
sub-micron  sized  protrusions  on  the  electrode  surface,  which 
act  to  "focus'*  the  electric  field  and  current  stream  lines. 
Primary  consideration  here  is  with  the  anode. 
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B.   REVIEW  of  PREVIOUS  SORK 
1 .   Physics 

Most  work  on  sheath  phenomena  is  embodied  in  "probe" 
theory  investigations,  where  the  affects  of  the  electrode 
(probe)  locally  disturb  the  quiescent  plasma.   Such  work  is 
relevant  to  this  study  since  the  anode  is  essentially  a 
heavily  biased  probe  in  contact  with  the  plasma,  which  is 
quiescent  within  the  sheath  region. 

Analyses  of  electrostatic  probes  in  collision  domi- 
nated plasmas  wece  originally  done  by  assuming  that  a  quasi- 
neutral  diffusion  controlled  region  extended  to  within  one 
nean  free  path  of  the  probe  where  it  was  matched  to  the  edge 
of  the  fr.ee  fall-sheath.   No  provision  was  made  for  a  tran- 
sition region  between  the  two  regimes  [Ref.  4]. 

Shultz  and  Brown  [Ref.  5]  studied  the  transition 
from  a  collisionless  to  collision  dominated  sheath  through 
charged  species  collection  with  a  spherical  probe.   Postu- 
lating different  physical  models  for  cases  ranging  from  a 
collisionless  sheath  to  a.  sheath  with  many  collisions, 
semi-empirical  formulas  were  obtained  in  each  instance. 
Several  other  authors  attempted  matching  the  solution  in  the 
collision  dominated  region  directly  to  that  of  the  free-fall 
sheath  using  variational  principles. 
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The  previously  mentioned  techniques  can  be  criticized, 
because  in  each  of  these  methods  something  was  forced  to 
fit.   In  other  words,  at  least  one  derivative  was  disconti- 
nuous at  some  artificially  imposed  boundary,  thereby  raising 
doubt  as  to  the  validity  of  the  results. 

The  first  systematic  analysis  of  probes  in  collision 
dominated  plasmas  were  carried  out  by  Cohen  [Hef.  6]  and  Su 
and  Lam  [3ef.  7],  who  assumed  the  continuum  equations  valid 
throughout  the  plasma,  including  to  the  probe  surface.   As  a 
consequence  of  these  analyses,  both  the  sheath  and  the  qua- 
sineutral  regions  appeared  naturally  from  the  diffusion 
equations  in  the  limit  as  the  ra-tio  of  probe  size  to  Debye 
length  approached  infinity.   The  continuum  model  removed  the 
concern  about  forcing  sets  of  equations  to  match  at  the  edge 
of  an  arbitrarily  assumed  sheath.   These  authors  restricted 
the  problem  to  an  isothermal  plasma,  but  as  will  be  shown 
here,  the  electron  temperature  may  vary  significantly  in  the 
neighborhood  of  the  probe  for  most  discharges  of  interest. 

Radbill  [Ref.  8]  guasilinearized  the  continuum  equa- 
tions, and  solved  with  numerical  integration,  obtaining 
results  that  agreed  with  those  of  Cohen  where  comparison 
could  be  made. 
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In   the  previous  analyses,    not   only    was   the    isother- 
mal   plasma   assumed,    but   the    plasma   was   so   slightly   ionized 
that   only   charge- neutral  interactions   needed   to   be   taken 
into   account.      However,    for    highly   nonequilibrium    plasmas, 
charge-charge  collisions   play  an  important   role  [Refs.    9  and 
10].      The  effects   of  collisions   between   charged   particles 
manifest  themselves   as  volume  ionization,    volume   recombina- 
tion  and  the  dependence   of    the    transport   coefficients  on 
charged   particle   density   and  temperature    and,    therefore   on 
position.      Barad   [Ref.    11]   cites  a   specific   example, 

"If   the  electrons   have   a    temperature  of    2500    °K,    then   the 
ionization  fraction   equal   to   about    10~s    will  result   in 
charge-charge   collisions    occurring    one-tenta   as   frequently 
as  charge-neutral  collisions.      It   is   shown...    that    for 
such    a  degree   of  charge-charge   interaction,   significant 
effects  are   felt." 

Chung,    Talbot   and   Touryan  [Ref.    1],   in   an    extensive 
review   of   probe   studies,    state   that   no   general   solution   is 
available   for  determining   charge  density    and  species   temper- 
ature  for    probes    small   in   comparison   to   the   hydrodynamic 
boundary   layer   thickness. 

Dolson   [Ref.    12]   investigated   the   nature   and   extent 
of   the    voltage   drops   in   the   vicinity   of  magnetohydrodynamic 
(M9D)    non-emitting  electrodes,    in   particular  the   losses 
attributable   to   the   sheath.      Under   his  conditions    the   nonex- 
istence  of   a  one-dimensional  sheath   solution  is   shown,    and   a 
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computer  model  with  two-dimensional,  periodic  active  sites 
representing  a  flat  plate  electrode  is  developed.   With  this 
model  the  effects  of  a  magnetic  field  and  Joule  heating  are 
studied  and  the  results  are  compared  with  experimental 
observations  reported  07  Argryoupolos  (1973) ,  and  Son ju  and 
Teno  (1974),  [Refs.  13,  14]. 

An  extension  of  Dolson* s  work  is  reported  in  1980 
[Hef.  15],  for  application  to  ionizar/sustainer  discharges. 
The  role  of  electron  pressure  is  introduced  along  with  net 
ionization  and  recombination,  which  proves  necessary  in 
establishing  truely  boundary-layer  behavior  for  the  charged 
species  profiles.   Additionally,  a  modified  Newton-Raphson 
computational  procedure,  which  enables  much  faster  solution 
convergence  is  formulated.   That  study  demonstrates  solu- 
tions with  non-active  wall  densities  depressed  from  the  free 
stream  density,  a  result  Dolson  did  not  accomplish.   Results 
were  for  10  volts  of  electrode  potential  fall,  an  increase 
from  5.0  volts  reported  by  Dolson. 
2.   Computational  Methods 

The  fundamental  difficulty  in  the  development  of 
computational  models  for  the  prediction  of  plasma  flow  phe- 
nomena lies  in  the  coupling  of  equations  governing  both  the 
gasdynamic  flow  and  the  electromagnetic  (EM)  fields.   Since 
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the  flow  fields  can  affect  the  EM  fields  (and  vice-versa), 
the  nonlinear  partial  differential  equations  must  be  solved 
simultaneously  and  the  formulation  of  the  problem  becomes 
complicated.   In  the  present  study  the  flow  properties  are 
decoupled  from  the  the  charged  species  conservation  equation 
and  Poisson's  Equation. 

The  approaches  used  in  the  solution  of  plasma  flow 
phenomena  are  quickly  reviewed  below  so  as  to  point  out  the 
background  of  the  the  approach  in  this  study.   This  problem 
is  typical  to  those  dealing  with  plasma-boundary  layer  phe- 
noaena.   Approaches  to  solving  the  general  plasma  flow  prob- 
lems have  consisted  of  the  following; 

(a)  One  dimensional  "influence  coefficients." 

(b)  One-dimensional  closed  form  analytic  solutions  to 
ordinary  differential  equations. 

(c)  Similarity  transformation  and  closed  form  solutions 
to  a  modified  differential  equation  in  one-  or  two- 
dimensions. 

(d)  One-,  two-,  or  three-dimensional  numerical  solution 
to  the  uncoupled  steady  plasma  flow  phenomena 
equations. 

(e)  One-,  two-,  or  three-dimensional  numerical  solution 
to  the  coupled  steady  plasma  flow  phenomena  equations 
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using  either  the  magnetic  induction  equation  or  the 
elliptic  equation  governing  the  electric  fields. 
(f)   One-,  two-,  or  three-dimensional  transient  numerical 
solution  to  the  coupled  plasma  flow  phenomena 
equations  using  either  the  magnetic  induction  or  the 
elliptic  equation. 
Methods  (a)  through  (c)  are  exact  but  are  restricted 
to  a  narrow  range  of  problems.   Approaches  (d)  and  (e)  pro- 
vide rapid  convergence  to  a  steady  state  solution  but  do  not 
consider  the  transient  behavior  between  the  gasdynamic  and 
electromagnetic  fields,  or  often  the  strong  coupling  between 
the  two  sets  of  fields.   Method  (f)  is  most  desirable  since 
it  treats  a  wide  range  of  problems,  but  suffers  from  large 
computational  times  [  Hef .  16].   The  approach  of  this  study 
is  characterized  as  a  method  (e)  approach. 

C.   OBJECTIVES  OP  THE  PRESENT  STUDY 

The  objective  of  this  study  is  to  demonstrate  a  numeri- 
cal procedure  which  can  model  the  physical  behavior  of  a 
conducting  fluid,  as  would  be  generated  in  the  after-glow  of 
an  electron  beam  discharge  laser,  under  the  influence  of 
electric  forces.   Calculations  are  performed  in  two  dimen- 
sions, with  a  flat  plate  non-emitting  electrode  (anode) 
which  has  periodic  active  sites.   These  active  sites 
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facilitate  current  constriction  required  to  satisfy  charge 
continuity,  Ohm*s  law,  and  Poisson's  equation.   The  physical 
model  is  idealized  in  Fig.  1. 

The  controlling  equations  are  solved  for  the  potential, 
the  charged  species  density  distributions,  electron  tempera- 
ture and  current  in  the  vicinity  of  the  electrode.   The 
resulting  computer  programs  generate  sheath  and  ambipolar 
regions  in  a  self -consistent  fashion,  using  the  same  set  of 
equations  throughout  the  field.   The  size  of  the  sheath  and 
the  voltage  drop  attributable  to  its  shielding  effects  are 
investigated. 

This  work  does  not  consider  those  phenomena  that  lead  to 
the  arc  discharge,  which  results  from  thermal  instabilities. 
The  constriction  at  the  anode  is  purely  that  required  to 
satisfy  continuity. 
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Two-Dimensional  Model  with  Periodic  Current 
Constriction. 
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II.       PROBLEM    FORHUIsATION 

A.       GOVERNING    EQUATIONS 

The    governing   equations    that   describe   the  electrode 
region    portion  of   the   plasma,   consist   of   Poisson's   equation, 
species   conservation   and  electron   energy    equation,    in   addi- 
tion  to    the   overall   continuity,    momentum   and  energy   equa- 
tions of  the   plasma.      These    overall   plasma   equations    are 
essentially   those   of   the   background   gas   since  the    plasma  is 
weakly    ionized    (less   than   0.01    %)  . 

We    have  [ Ref .    1 ] ; 


Poisson  equation... 
electron  conservation... 


O) 


IX- Vr^  -  V-[pe  (%?»*  ~  fr/cW>]  *     ^  (2) 


ion   conservation... 


ll-7A;-7-[D;(||7pi  +  ^r.^^)l  •     h< 


(3) 
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In  gas  discharges,  the  electron  energy  distribution  is 
known  to  be  non-Maxwellian.   In  order  to  allow  for  non-Max- 
wellian  distributions  and  still  retain  use  of  a  temperature 
concept,  an  empirical  fit  to  the  calculated  relation  of 
electron  temperature  as  a  function  of  E/n  is  used  in  the 
representation  of  the  electron  temperature  (see  Sec.  III. A) . 
The  concept  of  a  "two-temperature"  plasma  is  preserved  but 
no  restrictions  are  made  on  the  energy  distribution  function 
of  the  electrons.   In  addition  to  the  assumption  of  a  func- 
tional dependence  of  the  electron  temperature  upon  E/n  only, 
there  is  no  representation  of  an  external  magnetic  field 
and,  in  particular  the  internal  or  self-generated  magnetic 
fields  are  not  considered.   The  overall  gas  temperature  is 
assumed  constant  and  in  equilibrium  with  the  plasma  ions. 
The  bulk  energy  equation  is  assumed  to  have  little  bearing 
on  this  study,  since  there  is  negligible  gas  heating  as  dis- 
cussed in  Section  III. A.   This  effectively  uncouples  overall 
energy  consideration  from  the  problem. 

Tire  electrodes  serve  the  purpose  of  sustaining  the  E/n 
to  pump  the  laser  medium.   No  account  is  made  of  an  external 
E-beam  or  other  source  of  primary  ionization.   Simply,  the 
medium  may  be  considered  to  exist  in  an  afterglow  state. 


22 


The    production   term    (i.e.,    the   right   hand  side)    in  Eqns. 
(2)    and    (3)     has   been   described    utilizing   measured   and  theor- 
etical   results   of  the  rate   coefficiants  for   Nitrogen   ioniza- 
tion and  two-,    and  three-body  recombination.      The    production 
term   is    expressed  as; 

where   "^;"    and   "  &\  ■    are  respectively    the    ionization   and 
recombination   rate  coefficients. 

3.       BOUNDARY    CONDITIONS 

In    the   customary    boundary- value    problem   in   fluid    mechan- 
ics,   a    two-dimensional   flat    plate   is   assumed.      A    similar 
two-dimensional  Cartesian   description   is   used  in    the  present 
work,    which    has   been    found    to  be  the    minimum  suitable 
description    for    the   problem   posed,    [Refs.    3,12,15].      No 
one-dimensional   solution   exists   for   the   frozen   flow  of 
charges    through  a  collisional  sheath,    since  a  Cartesian  geo- 
aetry  in  one   dimension  is   devoid  of   the  necessary    geometric 
decrease  of   current    density   away   from   the   electrode,    as 
there   would   be   in   a   cylindrical   or   spherical  geometry. 

Convection    is    considered  negligible,    as   can    be   demons- 
trated  by   comparing   flu-id   velocity      with   the  drift   velocity. 
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Electron  and  ion  densities  are  determined  at  the  plasma 
edge,  knowing  tha  state  of  the  medium.   A  current  balance 
can  be  used  to  establish  the  species  densities  resident  at 
the  active  portion  of  the  anode.   An  electrode  voltage  drop 
is  determined  from  the  properties  of  the  medium;  e.g.,  the 
effective  ionization  energies  of  the  gas  in  particular. 
Since  current  constricts,  there  must  be  regions  of  the  flat 
plate  that  receive  no  current.   The  inactive  portion  of  the 
flat  plate  merits  attention,  especially  in  the  regions  adja- 
cent to  the  active  anode  site.   Species  concentration  along 
the  inactive  portion  of  the  wall  can  vary,  depending  on  the 
reactivity  (catalytic  properties)  of  the  wall.   There  is  a 
balance  that  can  be  established  between  the  voltage  varia- 
tion normal  to  the  wall,  and  the  species  concentration  at 
the  wall  together  with  the  species  variation  normal  to  the 
wall.   A  zero  concentration  of  electrons  and/or  ions  at  the 
wall  is  an  acceptable  boundary-condition;  which  is  to  say 
that  the  wall  acts  as  a  perfect  sink  for  the  charged  spe- 
cies.  A  detailed  kinetic-theory  analysis  shows  that  there 
is  in  fact  a  lower  bound  to  the  electron  and  ion  concentra- 
tions at  the  wall  [ Ref .  17].   In  any  case,  there  is  no  cur- 
rent into  or  out  of  the  inactive  wall,  (see  Fig.  1)  .   This 
is  discussed  later  in  Sect.  IV. 
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The  remaining  portion  of  the  boundary  conditions  refer 
to  the  free  stream  medium,  and  are  easily  established. 

C.   PARAMETERS  OF  THE  DISCHARGE  MEDIUM 

The  application  of  this  study  is  for  molecular  Nitrogen 
at  one  amagat#  where  the  gas  temperature  is  273  °K.   Other 
pertinent  parameters  follow; 

Electric  Pield:  1.5  to  15  KV/cm. 

21ectron  Density:  10«-»7  to  10«-*«  m~3. 

Electrode  Voltage  Drop:  35  volts. 
The  electron  energy  relation  is,  as  mentioned,  a  curve  fit 
from  experimental  results  [ Ref .  18],  carried  out  by  R.W. 
Compton  and  D.J.  Sutton  (1952).   Since  cross-section  data 
for  electrons  in  Nitrogen  discharges  are  available  [Ref. 
19  ],  the  electron  temperature  and  diffusion  coefficient  are 
readily  determined,  and  representable  as  a  function  of  E/n. 

The  ion  temperature  is  assumed  to  follow  that  of  the 
neutral  gas.   An  electrode  potential  fall  of  about  35  volts 
may  be  assumed  since  it  is  the  effective  ionization  energy 
for  Nitrogen  gas  and,  therefore,  the  anode  fall  [Ref.  20]. 
The  cathode  is  more  complicated  to  describe  because  of 
the  requirement  of  surface  electron  emission.   The  cathode 
fall  in  cold  cathodes  is  of  the  order  of  250  volts,  and  is 
considerably  less  in  thermionic  cathodes.   A  proper 

25 


III.  PROBLEM  SOLUTION 

A.   PROBLEM  SIMPLIFICATION  (PHYSICS) 

There  are  a  number  of  simplifications  that  are  incorpo- 
rated in  the  problem  description.   Namely; 

•  Steady  state  conditions  (within  the  pulse  time). 

•  No  magnetic  fields. 

•  Negligible  convection. 

•  Negligible  Joule  heating,  Tt-  =  T6=  constant. 
However,  specific  account  is  made  of; 

•  Ionization/recombina tion. 

•  Electron  pressure. 

These  are  discussed  in  part  in  this  section. 
1 .   Convection 

The  ambipolar  region  is  the  transition  region  from 
the  sheath  to  the  undisturbed  plasma  and  can  perhaps  span 
the  boundary  layer.   This  fact  requires  some  consideration 
of  the  presence  of  convective  effects. 

From  the  ion  species  conservation  eguatio.n  (Eq.  3)  , 
for  example,  we  know  that  in  the  ambipolar  region  the  last 
term  of  the  left  hand  side  is  small  compared  to  electric 
conduction.   Now  we  need  only  to  compare  the  convection  with 
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the  conduction,  or  simply  the  fluid  velocity  with  the  drift 
velocity. 

v  =  y^R   «  0.1  to  1.0  m/sec 
where  the  Reynolds  number  (Re)  is  10*  to  10*  and  the  free 
stream  velocity  is  10  to  100  m/s.   Now  the  drift  velocity 
may  be  approximated  by 

vd  «»  JA^-   10  to  100  m/s 
for  E  -  105  to  10*  V/m  and  the  mobility  /d  =  10-* 
az/s/volt.   Clearly,  in  the  presence  of  a  relatively  strong 
interelectrode  field,  the  contribution  of  ion  convection  due 
to  a  cross  flow  should  be  negligible.   The  geometrical 
orientation  of  the  discharge  with  respect  to  the  primary 
flow  can  thus  become  an  important  consideration. 
2.   Density  Sounder.!  Layer  Formation  Time 

The  charged-species  density  profiles  change  in  a 
more  complicated  manner  than  either  the  voltage  potential  or 
the  electric  field.   The  magnitude  of  the  ion  and  electron 
densities  may  change  appreciably  within  the  sheath;  more- 
over, fractional  analysis  is  risky  because  inflections  are 
present  in  the  profiles.  k    discussion  pertaining  to  the 
stability  of  the  ambipolar  region  and  its  boundary  layer 
nature  is  presented  in  detail  by  piblarz,  et  al.  ,  [Ref.  15]. 
However,  Appendix  &    shows  that  the  characteristic  length  of 
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the  ambipoiar  region  is  of  the  order  of  the  characteristic 
sheath  length.   Letting  L=5x10~5  meters  and  the  drift  veloc- 
ity as  before,  vd  =  10  to  100  m/s,  a  characteristic  time  for 
the  formation  of  the  density  boudary  layer  would  be  0.5  to 
5.0  microseconds.   It  is  reasonable  that  the  ambipoiar  dif- 
fusion to  the  wall  can  establish  itself  well  within  the  50 
microsecond  pulse  period  of  the  electrodes. 
3-   Significance  Of  Joule  Beating 

For  the  present  problem,  joule  heating  can  be  shown 
not  to  be  significant.  In  the  vicinity  of  the  active  site, 
the  average  energy  deposition  rate  per  unit  volume  is  shown 
to  be  approximately  Ux108  watts/m3.  The  temperature  of  the 
Nitrogen  gas  at  one  amagat  would  bs  raised  approximately  1°K 
in  2.6x10-*  sec.  allowing  a  50x10-*  sec  pulse  period  would 
raise  the  temperature  of  that  region  by  approximately  20°K, 
or  7%. 

Assume  that  the  heat  generated  at  the  anode  spot  is 
carried  away  by  thermal  diffussion  (conduction).   Then  for 
Nitrogen  gas  (Jc=0.0267  W/m-°K,  /  =1.138  kg/m3,  and  cp=1043 
J/kg-°K)  the  thermal  diffusivity  is; 

Ots  -pr  =  2.25x10-3  m2/sec 
A  typical  diffussion  time  would  be  (for,  A3  =4.4x10-5  m)  ; 

?de\/o(z,   0.85  millisec 


'.J 
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The  "flush  time"    is    the   time    required   for   clearing 
the   lasing   cavity   of   residue   gases   in    preparation    for   the 
next   pulse,    and   can   be   estimated   as   the   gas  flow   rate 
divided    by   the   flow   length.      Assuming   v   =    100   m/sec   and   L  = 
10    cm,    the    flush   time   would    be; 

FT   =   0.1/100   =   1    millisec 
Other  factors  considered,    the  inter-pulse   period    (time  bet- 
ween  pulses)    would  be   larger,    but   even  in   a   one   millisecond 
period,    there   is    just   enough  time   to   cool    the  anode   spot. 

The   analysis   so  far   has  assumed  a    flow   velocity   at 
the   surface  to   be  the  same    as   the   bulk   flow   in    the   laser 
cavity.       However    a   boundary    layer   exists —    probably   turbu- 
lent,   with    a   laminar   sub-layer   adjacent   to    the   electrode 
for; 

Rc    =    XV^/,?    =    (0.1)  (100)/(  15.53x10-*)     =   6.44x105 
The   velocity  boundary   layer   thickness    would   be   2.55    millime- 
ters.     Assuming   a    1/7-th   power    profile,    the   velocity    at  a 

sheath   length  from   the  electrode   would   be: 

i/i 
v  =    V(Aj/$)     -    100(4. 4x10-s/2. 55x10-3)  V7   =    55    ffi/Sec 

A  characteristic  local  flush  time  would  be  approximately  2 
milliseconds,  and  much  higher  for  a  laminar  velocity  pro- 
file.  These  flush  times  can  be  quite  large  as  the  require- 
ment to  cool  the  Joule  heated  gas  near  the  anode  spot  is 

30 


increased.   The  design  flush  time  based  on  other  factors  may 
infact  be  insufficient  to  clear  enough  of  the  Joale  heated 
gas  from  the  vicinity  of  the  anods  spot,  enabling  a  build  up 
of  heated  gas,  which  could  lead  to  local  enhanced  ionization 
of  the  gas. 

4.   Electron  Temperature 

The  electron  temperature  (Te )  is  assumed  to  be  a 
linear  function  of  E/n  as  presented  in  Fig.  2,  for  the 
region  of  interest,  namely  E/n  =  0.5  to  5.0x10~2°  V-m2.   For 
the  range  of  2/n  considered,  it  is  adequate  to  take  the 
electron  temperature  as; 

Te  =  [A«log(E/n)  ♦  B]«T0  (5) 

where; 

A  =  12.  1 
B  =  38.5 

E/n   is   normalized    to    10+20    V-m2 
The   diffusion    coefficients   were  determined   from 
Brown  [Ref.    18].      The   drift    velocities   were   obtained    from 
plotted   results,    and   since  jU^-Vj  /E,    and   D5  =;Li,kT5/e,    (Ein- 
stein's  relation),    the  diffusion   coefficients  are; 

Dg    =   0. 0665    m2/sec 
D.    =   8.0x10-*   m2/sec 
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B.       FORKING    EQUATIONS 

Since  convection   is   negligible,    and   assuming   that   the 
sustainer   operates   in  the   afterglow  of   an   electron    beam,    the 
governing  equations,    namely    Eqs.    1,2    and   3,    can   be   written 
as    follows; 

V24  =    c,  (ne-  n;)  (6) 

-V»[Vne+  n^V^/d   -   n^W/d]  =C2.ne  (7) 

-7»(7n;  ♦  £n(-V0]  =    c,nd  (8) 

where; 

c,=    e^/kT^  ct=  As2n«^/De  c3=   c2>D4/Dt- 

The    equations   may   now    be  applied   to   flat   plata    formula- 
tion  as    depicted    in   Fig.    3.      It   is   convenient  to   drop   the 
ivmi    for    simplicity,    but  it    is  to  be    understood   that   all 
variables  have   been   suitably  non-dimensionalized. 

We    have   Poisso^s   equation; 

ft  *  S^ '  *iH'mi  <9' 
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electron  conservation  equation; 


ion  conservation; 


c.K<MO-*[#g*$£Hfe+&]-   c^       (11> 


C.   NUMERICAL  METHODOLOGY 

The  nonlinear  terms  creat  problems  in  the  computer  anal- 
ysis of  the  electron  and  ion  species  conservation  equations, 
(Sqs.  10,11).   The  Jacobi  method  includes  all  nonlinear 
terms  on  the  "right  hand  side";  i.e.,  external  to  the  coef- 
ficient matrix.   Convergence  to  a  solution  is  possible  if 
these  non-linear  terms  change  slowly  over  each  iteration. 
It  was  found  by  Dolson  [Ref .  12  ]  that  the  Jacobi  method  was 
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in   fact    unstable    for   the   present  set   of   equations   and 
conditions.      As   a   consequence  of   the   failure  of   the   Jacobi 
method,    a  quasi-Jacobi  method  was   implemented,    in    which   an 
estimate  for   each   of   the   solution  variables    {4> ,   nfe ,    v.L)    was 
computed.      Hhen   the   product   of   two   variables   is   encountered, 
one   variable   is   treated  as   a  constant   coefficient    for   each 
iteration.      This    means  that    the   non-linear   terms   are 
retained  in   the  coefficient    matrix.      The   "constant"    coeffi- 
cients   are   updated  after   every   iteration,    thus  altering  the 
coefficient   matrix.      The   conventional   Jacobi   method   was 
found  to   converge   only   for    low   values   of    voltage   potential 
at   the    electrode,    whereas    the  quasi-Jacobi    procedure    pro- 
vided  converged   solutions   for  approximately   5  volts.      A   dis- 
cussion  of    these    methods   is    provided   in   Appendix  C    as   back- 
ground   for   the  development    of  the   modified   flewton-Raphson 
(MMR)    method,    which   is  used    in   the   solution   to   Eqs.    1 
through    4. 

The   nonlinear   coupled   partial   differential   equations 
(PD2) ,    Eqs.    9,10    and    11    are    linearized   to    provids   a   set  of 
linear    algebraic   equations,    which   must   hold   throughout   the 
domain    described   as   follows.      The  domain   as  illustrated   in 
Pig.    3,    is   divided  up   into    equal   sided  sub-areas   whose   ver- 
tices  are  the  computational    nodes  at    which   the   linearized 
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algebraic  equations  must  hold.   The  wall  is  represented  by 
seven  (7)  nodes,  with  the  leftmost  node  being  the  location 
of  the  active  site.   The  "y"  direction  extends  from  this 
active  site  along  the  wall.   The  »x"  dimension  perpendicular 
to  the  wall  is  represented  by  25  nodes  starting  from  each  of 
the  wall  nodes.   The  solution  to  the  PDE^  is  represented  by 
the  vector  z(x,y)  ,  or  on  the  computational  mesh  as  z (i, j) . 
The  boundary  conditions  have  been  discussed  previously 
and  are  presented  in  Pig.  3.   The  boundary  condition  that 
there  is  no  net  current  into  the  nonactive  wall  {J*  =0)  is 
satisfied  by  the  electron  contribution  to  current  only. 
This  is  justified  since  it  can  be  shown  that  the  electron 
contribution  is  approximately  three  orders  of  magnitude 
greater  than  the  ion  contribution  to  the  total  current. 
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Boundary  Conditions  for  Computational  Method. 
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IV.  RESULTS  AND  CONCLUSIONS 

A.   GENERAL  DISCUSSION 

Computer  solutions  are  reported  here  for  three  cases 
presented  in  Table  I.   The  active  site  voltage  was  35  volts 
in  each  of  the  three  cases.   Included  in  this  section  are 

TABLE  I:  Computer  Solution  Conditions 


CASE  ne(wall)  n^  ( wall)  ne  (anode)   n£  (anode) 

I  1.0  0.0  J-matched    0.0 

II  "float"  0.0  not  J-matched   0.0 

III  "float"  "float"  not  J-matched  0.0 


also  some  pertinent  results  obtained  at  10  volts  to  illus- 
trate the  effects  the  net  production  term  has  on  the  conver- 
gence to  a  boundary  layer  behavior  for  the  species  density 
and  free  stream  electric  field. 

The  solutions  are  depicted  as  oblique  presentations  of 
the  potential,  species  density,  space  charge,  net  produc- 
tion, electric  field  and  Joule  heating  distributions  over 
the  "x-y"  coordinates  defined  in  the  problem  formulation. 
Also  included  are  two-dimensional  plots  of  various  terms 
along  a  "cut"  or  line  perpendicular  to  the  electrode  wall, 
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starting  at  the  active  site,  with  some  starting  at  the  wall 
along  the  inter-nodal  line  of  symmetry  (see  Fig.  3)  .  The 
solutions  are  organized  as  Case  I, II  and  III  results,  will 
be  presented  as  Pigs.  12,13  and  14  respectively  at  the  end 
of  this  section.  The  discussion  of  these  results  proceeds 
with  additional  figures  which  are  interspersed  in  the  text 
appropriately. 

Table  I  presents  various  electron  and  ion  boundary  con- 
ditions along  the  wall  some  of  which  are  "floated."   By  this 
teem  is  meant  that  upon  the  completion  of  each  numerical 
iteration,  the  electron  (and/or  ion)  density  at  the  nonac- 
tive  portion  of  the  wall  is  re-established  by  assigning  it 
the  value  of  a  near  neighbor  away  from  the  wall.   This  pro- 
cess is  illustrated  in  Fig.  4. 


t 


V  I- 


-•x 


INITIAL  CONDITIONS 


I        V  *-? 


AFTER  ITERATION 


V  l 


'FLOATED' 


Figure  4 
Description  of  "floating  the  Wall  Boundary  Condition 
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The  usual  elliptic  boundary  value  problem  has  well 
possed  boundary  conditions  from  a  priori  information.   How- 
ever, the  nature  of  the  problem  in  this  study  includes  very 
little  information  as  to  the  correct  species  wall  density. 
It  seems  reasonable  to  set  these  densities  to  zero  along  the 
inactive  portion  of  the  electrode  as  previously  mentioned  in 
Sec.  II. B.   In  this  study  the  wall  boundary  conditions  are 
"lossened"  to  facilitate  "naturally  occuring"  and  unpresum- 
ing  species  density  profiles.   The  nature  of  the  problem  in 
effect,  produces  a  self-consistent  boundary  condition. 

B.   EFFECTS  OP  NET  PRODUCTION 

Earlier  work  [Ref.  15]  reported  limited  results  for  a  10 
volt  active  node  site,   solutions  for  cases  with  and  without 
the  inclusion  of  the  net  production  term  have  profound 
influences  on  the  appearance  of  the  charged  species  density 
profiles.   Figures  5,6  and  7  illustrate  the  before  and  after 
effects  of  including  the  net  production  term.   In  both  cases 
the  wall  boundary  values  for  the  species  density  was  ini- 
tialized at  the  free-stream  value  (i.e.,  1.0),  then  allowed 
to  "float'1  as  described  previously.   The  addition  of  the  net 
production  term  brings  about  the  boundary  layer  behavior  of 
the  densities  of  the  electrons  and  ions  is  illustrated  in 
Fig.  7.   Efforts  to  attain  converged  solutions  with  the  net 
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production    tern    set    to   zero    proved    fruitless   for    active   site 
potentials    exceeding    10   volts. 


Electron    Density 
(a) 


Ion    Density 
(b) 


Figure   5 


(a)    2lectron   and    (b)    Ion  density    perspectives    with 
the  net   production   term   set   to    zero.     (10    volts) 
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Pigure  6 

(a)    Electron   and    (b>    Ion  density    perspectives    with 
the   net   p reluct ion    term  coupled.     (10    volts) 
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Figure   7 

Influence   of    the   net   production    term   on    the 
density    profiles    along   a    line    perpendicular   to    the    wall, 
extending   from    the   active   node    site,     (10    volts). 
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C.   DISCUSSION  OP  RESULTS 
1.   Voltage  Potential 

The  voltage  potential  of  the  active  node  site  ("the 
anode")  in  each  of  the  cases  presented  was' 35  volts.   Fig- 
ares  12-a,  13-a  and  14-a  will  clearly  illustrate  sharp  vol- 
tage spikes  at  the  anode,  and  the  "Laplacian"  nature  of  the 
field  far  from  the  anode.   The  potential  along  the  inactive 
wall  varied  from  12  volts  at  the  farthest  position  from  the 
anode  to  17  volts  just  next  to  the  anode.   Similar  profiles 
of  the  electric  potential  along  the  wall  for  the  other  cases 
are  presented  in  Fig.  8.   Figures  12-i#  13-i  and  14-i  will 
show  the  voltage  potential  variation  along  a  line  normal  to 
the  wall  extending  from  the  anode  site. 

Figure  9  illustrates  a  two-dimensional  cut  of  the 
voltage  potential  along  the  inter-nodal  symmetry  line.   Here 
it  can  be  seen  that  the  variation  from  the  wall  is  shallow. 
It  is  important  to  point  oat  that  the  electron  temperature 
is  directly  related  to  the  gradient  of  the  voltage  potential 
through  Eg.  5.   This  fact  has  certain  consequences  which  are 
discussed  in  the  next  section. 
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Pigure  8 

Voltage  potential  along  the  wall  boundary  (electrode) 
far  CASES  I,  II,  and  III. 
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Figure    9 

Voltage   potential  variation   normal    to    the    wall, 
along   the    inter-nodal  symmetry    line,     (CASE    I) 
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2.   Electron  Enerqv. 

The   electron   energy    will  be   presented  in   series   -b 
of   Figs.    12,13   and    14 ,    as   oblique  perspectives.      The   anode 
site  temperature   for  Cases    I,   II  and   III    are   2.1,    2.1   and 
1.9    e7   respectively.      The   description   of   the  electron   temp- 
erature   is   that   of   a   linear    function   of   the   logarithm   of  the 
electric  field.      Any    portion  of   the    potential   field   that   is 
"level"    (grad  <f>  =0)    will   necessarily    generate  a   low   or  even 
negative    (see  Fig.    14-b)    electron  temperature.      according  to 
Fig.    2   and   the   discussion   in   Sec  II. A. 4,    a   negative   tempera- 
ture   would   occur    below   an   E/n  of  0.04x10-2°   7-cm2.      This 
problem    can    be  removed  by   setting   a    minimum    value   to   the 
non-dimensional  temperature    field,    nominally  0.13,    corres- 
ponding   to    the   minimum  at   the  low  end   of    E/n   in   Fig.    2. 

The   parameters  of   the  discharge    (Sec.    I.C)    require 
an   electric   field   of   2.7x10s   V/m    (or,    log£ E]=5. 43)    at   the 
free   stream   boundary,    which    corresponds  to   an  electron   temp- 
erature  of    0.905    eV    (<2=1.)  .      The   temperature  fields    in  all 
cases  presented   do  not  quite  meet   the   free   stream   boundary 
requirements.      The  electric    fields   are   insufficient   to   ren- 
der  proper    results  to   the   electron   temperature   profiles. 
Thus   the   free  stream  temperature  boundary   requirements   are 
"forced"    by   slightly    inflating   the   influence  of   the   electric 
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field.   The  effect  of  forcing  the  temperature  field  in  this 
way  can  also  be  looked  at  as  "sliding"  the  linear 
representation  of  the  temperature  curve  to  higher  values  at 
constant  B/nr  (see  Fig.  2).   Figure  10  illustrates  a  case 
where  a  low  electric  field,  here  0.8  lcV/cm,  results  in  a  low 
electron  temperature  projecting  toward  the  free  stream 
(0=0.62).   The  low  temperature  is  adjusted  to  meet  the 
desired  free  stream  value  (0=1.0)  by  multiplying  the  elec- 
tric field  by  a  factor  of  3. 40.   Figure  10  also  illustrates 
an  alternate  nethod  of  meeting  the  free  stream  conditions  by 
truncating  the  variation  of  temperature  with  a  minimum 
value,  namely  the  free  stream  boundary  condition. 

Both  of  the  "forcing"  techniques  are  used,  but  for 
different  reasons.   As  long  as  the  adjustment  factor  does 
not  exceed  an  order  of  magnitude  change  in  the  electric 
field,  then  it  is  not  objectionable  in  light  of  the  simpli- 
fied representation  of  the  electron  temperature.   The 
adjustment  requirements  of  the  three  cases  reported  in  this 
study  are  between  3.0  and  5.0,  with  Case  I  using  a  factor  of 
10  to  see  the  effects  on  the  Joule  heating.   In  Case  I  the 
non-dimensional  temperature  is  artificially  high  at  1.4  eV 
at  the  free  stream. 
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Figure    10 

Electron   tenperature    variation    normal    to    the    wall,    along   the 
nodal   symmetry   line;    illustrating   two    methods  of    adjustment 
to    '•force"    the    free   stream  temperature    conditions. 
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The  electric  field  near  the  inter-nodal  symmetry  points 
of  the  wall  are  characteristically  2  orders  of  magnitude 
less  than  that  of  the  free  stream  conditions.   The  resulting 
very  low,  even  negative  temperatures  in  Case  III,  Fig.  14-b, 
is  "patched"  using  the  truncation  method  previosly  dis- 
cussed,  whether  these  low  temperatures  and  the  associated 
steep  derivatives  are  physically  reasonable,  merits  further 
investigation. 

The  temperature  fields  generated  and  presented  as 
Figs.  12r  13,  14  series  -b  are  the  final  result  depictions. 
However,  the  temperature  fields  used  in  the  solution  process 
were  adjusted  and  truncated  as  previously  discussed  to  main- 
tain solution  convergence.   An  example  of  a  "corrected" 
temperature  field  is  presented  in  Fig.  14-c. 

3.   ^lectris.  £ieii  And  Joule.  Hi.  at  in  a 

The  electric  field  and  Joule  heating  distributions 
are  generated  from  the  primary  voltage  and  electron  density 
solutions,  much  in  the  manner  as  the  temperature  distribu- 
tion was  determined.   In  fact,  the  Electric  field  and  the 
Joule  heating  distributions  are  similar  in  appearance  to  the 
temperature  distributions.   These  distributions  play  no 
direct  part  in  the  solution  scheme  and  thus  being  ancillary 
in  nature  were  not  "forced"  to  meet  any  a  priori  conditions. 
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The  Joule   heating    term   can  be    written  as; 

and  thus   is   a  rather   straight-forward   function   of     <t> ,   ne, 
and   n,- .      The   ion   contribution   is  shown   to   be  negligible, 
since   De  >>    Dc-  .      The   Joule   heat    at   the   anode   site   is   of   order 
10s   and   is    100   times    that   of  the  fres    stream   values.      Joule 
heating   effects   are   negligible    in   the   overall   problem,   as 
discussed  in   Sec.    II. A. 3,    obviating   the  requirement   for   an 
energy   eguation. 

4.   Species  Density  Distribution 

The  solutions  for  the  charged  species  density  dis- 
tributions are  presented  as  oblique  perspectives  in  Figs. 
12,13  and  14,  series  -e  (electron  density);  series  -f  (ion 
density);  and  series  -g  (space  charge  density).   Addition- 
ally the  variation  of  these  distributions  normal  to  the  wall 
along  the  active  node  symmetry  line  are  presented  in  the 
series  -j  and  -k  -figures.   The  oblique  perspectives  of  the 
densities  are  reversed  foe  clarity. 

The  electron  and  ion  densities  at  the  nonactive  wall 
are  a  subject  of  interest  in  the  literature,  to  which  little 
information  has  been  gathered  on  a  firm  basis.   The  results 
of  this  study  suggest  a  slightly  positive  space  charge  at 
the  inactive  wall.   Scrutiny  of  Case  III  results  reveals,  a 
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2.2x1016  m-3  positive  space  charge  number  density,  (or  3.5 
milli- coulombs) .   Recall  that  3ase  III  involved  "floating" 
both  of  the  charged  species  densities  at  the  wall.   In  fact, 
close  inspection  of  Case  I  and  II  reveals  a  positive  space 
of  similar  magnitude  but  displaced  slightly  from  the  wall. 
The  wall  in  these  cases  was  set  with  a  negative  space  charge 
as  the  boundary  condition.   The  positive  space  charge  along 
the  wall  is  most  predominant  in  the  vicinity  of  the  inter- 
nodal  symmetry  line,  and  whether  the  low  temperature  values 
in  this  region  have  significance,  remains  to  be  determined. 

The  active  node  site  (anode)  electron  density  in 
Case  I  (2.4x1017  a-3)  was  determined  so  as  to  match  the  free 
stream  current  to  the  current  constricting  at  the  anode.   In 
the  results  of  Cases  II  and  III  the  anode  electron  densities 
are  chosen  so  as  to  provide  a  smooth  variation  extending 
outward  from  the  active  site. 

5.   Ionization/Recombination 

The  net  ionization  term  is  utilized  in  all  cases; 
and  was  found  necessary  for  (1)  establishment  of  the  charged 
particle  density  distributions,  and  (2)  solution  convergence 
for  sheath  lengths  spanning  more  than  one  computational  node 
spacing.  The  net  production  term  is  depicted  in  Figs.  12,13 
and  14  series  -h.   There  is  net  ionization  throughout  the 
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solution  domain.   However  as  Case  III  demonstrates,  there  is 
a  factor  of  5  orders  of  magnitude  greater  ionization   within 
the  sheath  region  surrounding  the  anode  site  than  at  the 
free-srream  edge.   Cases  I  and  II  show  ionization  along  the 
nonactive  wall,  simply  because  the  ion  density  is  zero,  eli- 
minating significant  recombination,  (see  Eg.  4)  . 
6»   Current  Description 

Figure  11  is  typical  of  the  constricting  nature  of 
the  current  stream  lines.   The  potential,  electron  density 
and  the  electron  temperature  are  used  to  determine  the  "x-y" 
components  of  the  current  density.   Using  a  "nearest  neigh- 
bor" influence  principal,  a  local  derivative  can  be  deter- 
mined; from  which  a  streamline  can  be  "marched-out"  in  small 
increments.   The  profiles  of  all  three  cases  are  similar, 
the  fact  being  that  the  current  matching  discussed  previ- 
ously will  not  affect  the  current  streamline  profiles.   Fig- 
ure 11  clearly  demonstrates  the  satisfaction  of  the  current 
constriction  requirement  necessary  to  permit  solutions  to 
the  system  of  equations. 

The  current  into  the  anode  site  in  Case  I  is  matched 
to  the  current  entering  at  the  free-stream  boundary.   This 
specifies  the  electron  density  at  the  anode  site.   However, 
in  Cases  II  and  III  the  current  at  the  anode  was  only  three 
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times   the   free-stream  currant.      The   current   matching    problem 
does  not   deserve    serious   complaint,    since   the   problem 
description    is   simplified    tc  two-dimensions  and   the   current 
constriction   is   only   required  to   satisfy   continuity. 
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^^^^^^^^ 
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NONACTIVE  WALL 
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Figure    11 

Two-dimensional    depiction   of   the   current   stream   line 
constriction   toward  the   active   node  site. 
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Figure  12-a 

Oblique  perspective  of 
Potential  Field  (CASE  I) 


Figure  12-b' 


Oblique  persi 
Temperature  fi< 


•  ective  of 
ild  (CASE  I) 


Figure    1  2-c 

Oblique    perspective   of 
log   Electric   Field     (CASE    I) 


Figure    12-d 

Oblique    perspective    of 
log    Joule    Heat    (CASE    I) 
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Figure    12-e 

Oblique    perspective   of 

Electron    Density    (CASE    I) 

(reversed    perspective) 


Figure   12-f 

Oblique    perspective  of 
Ion   Density    (CASE    I) 
(reversed   perspective) 


Figure    12-g 

Oblique    perspective   of 
Space    Charge    (CASE    I) 


Figure    12-h 

Dbligue    perspective   of 
Met    Production    (CASE   I) 
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Figure  12-i 

Presentation   of    the  Voltage    Potential 

(CASE    I)    along  a    line  normal    to   the  wall, 

extending   from    the  active    node    site. 
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Figure    12-j 

Prasentation    of   Species    Density 

(CASE    I)     along  a    line    normal    to    the    wall, 

extending    from    the    active    node    site. 
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Figure    12-k 

Presentation    of   the   negative   Space    Charge    Density 

(CASE    I)     along  a    line    normal    to    the   wall, 

extending   from    the   active    node   site. 
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Figure    12-1 

Presentation    of   the   Electric   Field 
(CASE    I)     along  a    line    normal    to    the    wall, 
exteading    from    the    active    node   site. 
(logarithm   scaled) 
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Figure    12-m 

Presentation    of   the  Joule   Heating 
(CASE    I)     along  a    line    normal    to    the    wall, 
extending    from    the    active    node   site. 
(logarithm    scaled) 


60 


CO 

LiJ 

a 

CO 


CJ 
LU 

ci- 
cn 


1.0 


0.8 


0.6 


0.4 


0.2 


0.0 


!  I  I 


J !_ 


J I 


JL 


±_J I 


0.00 


0.05 


0.  10 


0.  15 


0.20 


DISTANCE     (MM 


exte 


Figure    12-n 

presentation   of    the    Species    Densities 
(CASS    I)    on   a   Line   normal   to   the   wall, 
nding   along   the    inter-nodal   symmetry   lin< 
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Figure    1  3- a 


Oblique    perspective  of 
Potential" Field    (CASE   II) 


Figure   13-fcf 

Dbligue    perspective   of 
Temperature   field    (CASE   II) 


Figure  1 3-c 

Oolique  perspective  of 
log  Electric  Field  (CASE  II) 


Figure    13-d 

Obligue    perspective   of 
log    Joule    Heat    (CASE    II) 
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Figure  1  3- e 

Oblique  perspective  of 

Electron  Density  (CASE  II) 

(reversed  perspective) 


Figure  13-f 

Oblique  perspective  of 
Ion  Density  (CASE  II) 
(reversed  perspective) 


Figure    1 3-g 

Oblique    oerspective   of 
Space    Charga    (CASE  II) 


Figure    13-h 

Dblique    perspective   of 
Net    Production     (CASE    II) 
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Figure    13-i 

Presentation  of    the    Voltage   Potential 
(CASE   II)    along   a    line    normal   to    the  wall, 
extending    from    the    active    node    site. 
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Figure    13— j 

presentation    of   Species    Density 
(CASE    II)    along   a    line    normal   to    the  wall, 
extending    from    the    active    node    site. 
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Figure    13-k 

Presentation    of   the   negative   Space   Charge   Density 
(CASS   II)    along  a    line    normal   to    the   wall, 

extending    from    the    active    node    site. 
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Figure    13-1 

Presentation   of  thfl   Electric    Field 
fCASF.   Ill    alona   a    line    normal   to    the  wall, 
^extending    from    the    active    node    site, 
(logarithm    scaled) 
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Figure    13-o 

Presentation   of  the  Joule   Heating 
(CASE   II)    along   a    line   normal   to    the  wall, 
extending   from    the   active    node   site, 
(logarithm    scaled) 
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Figure    13-n 

Presentation   of    the   Species   Densities 
(CASE    II)     on  a    line  normal    to   the    wall.. 
Indinq    along   the    inter-nodal  symmetry   lin« 
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Figure    H-a 
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Pigure   14-b 

of  Oblique    perspective   of 

III)  Temperature    field    (CASE   III) 


Figure    14-c 

Oblique    perspective   of   adjusted 
Temperature   Field    (CASE    III) 


Pigure  14-d 

oblique  perspective  of 
Joule  Heat  (CASE  III) 
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Figure  14-e 


Oblique  perspective 
Electron  Density 


of 
(CASE    III) 

(reversed   perspective) 


Figure    14-f 

Obligue  perspective  of 
Ion  Density  (CASE  III) 
(reversed   perspective) 


Figure    14-g 

Oblique    perspective   of 
Space    Charge    (CASE    III) 


Figure    14-h 

Oblique    perspective   of 
Net    Production    (CASE   IIJ) 
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Figure    1U-i 

Presentation  of    the    Voltage   Potential 

(CASE   III)     along    a    line   norm    to    the   wall, 

extending   from    the   active    node   site. 
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Figure    1U-j 

Presentation    of  Species    Density 

(CASE   III)    along  a    line    normal    to    the    wall, 

extending    from    the    active    node   site. 
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Figure    1U-lc 

Presentation    of   the   negative   Space    Charge    Density 

(CASE   III)    along  a    line    normal   to    the   wall, 

extending    from    the    active    node   site. 
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Figure    14-1 

presentation   of  the   Electric   Field 
(CASE   III)    along   a    line    normal   to    the   wall, 
extending    from    the   active    node   site, 
(logarithm   scaled) 
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Figure    14-i 

Presentation   of   the  Joule   Heating 
(CASE   III)    along   a    line    normal   to    the   wall, 
extending    from    the   active    node    site. 
(logarithm   scaled) 
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Figure    14-n 

presentation  of    the    Species   Densities 
(CASE    III)    on    a    line    normal    to    the    wall, 
extending    along    the    inter-nodal   symmetry    line 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 

A.   REVIEW  OF  ASSUMPTIONS 

In  the  development  of  this  work  several  assumptions  were 
made  which  the  solution  results  later  verified.   The  assump- 
tion of  negligible  convection  effects  is  straight  forward. 
The  charged  particle  density  boundary  layer  evolution  time 
turns  out  to  be  small  in  relation  to  the  pulse  period;  thus 
enabling  consideration  of  the  problem  as  a  steady-state  aft- 
er-glow problem.   The  fact  that  the  Joule  heating  of  the  gas 
is  shown  negligible,  frees  the  formulation  from  an  elaborate 
energy  description.   Since  the  electron  temperature  is  a 
simple  function  of  E/n,  the  problem  can  be  further  simpli- 
fied.  The  solution  proceeds  with  nominal  values  of  the  par- 
ameters of  the  discharge  as  discussed  in  Sec.  II.  B. 
Finally,  the  importance  placed  on  the  net  production  term  is 
verified  by  the  profound  changes  in  the  results  including 
the  fact  that  "ne"  is  necessary  for  the  establishment  of  the 
charged-species  boundary- layer  behaviour.   The  analysis  of 
App.  3  shows  that  an  assumption  of  "frozen  flow"  would  in 
fact  be  premature. 
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B.   PHYSICAL  CONCLUSIONS 

The  present  model  of  the  sheath  and  ambipolar  region 
evolves  from  the  assumptions  of  steady  state,  uniform  ion 
temperature  and  constant  diffusion  coefficients.   A  summary 
of  important  results  is  presented  in  Table  II. 

TABLE  II:  SUMMARY  OF  IMPORTANT  RESULTS 


SUMMARIZED 

RESULTS 

CASE 

I 

CASE 
II 

CASE 
III 

(units) 

sheath  length 

3.3 

3.5 

3.7 

10-s  a 

ambipolar  length 

4.5 

3.5 

3.7 

10-5   jg 

E^(f ree- stream) 

7.9 

3.9 

6.3 

10*  V/m 

Z,   (anode) 

1.9 

2.5 

2.2 

10*  V/m 

The  results  show  that  a  sheath  develops  as  anticipated 
from  theory  and  the  proposed  assumptions.   The  sheath  char- 
acteristic length  is  consistent  with  that  predicted  by  Debye 
shielding.   Debye  shielding  of  a  flat  plate  predicts  that 
the  sheath  length  M/\|"  is  U.4x10-5  meters.   For  a  point  (or 
spherical)  shielding  the  sheath  length  is  3.1x10-5  meters. 
The  present  study  involves  the  shielding  of  a  point  voltage 
source  situated  on  a  surface  charged  plate.   Since  both  two- 
and  three-dimensional  effects  are  inherent  in  the  formula- 
tion, the  resulting  sheath  lengths  are  quite  reasonable. 
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As   Table   II   shows,    the  size    of   the   ambipolar  region   is   of 
the   same  order   as   the   sheath  size.      Here   the  ambipolar 
length    is   defined  as   the   length  extending   from   the   edge  of 
the   sheath   to   the   point   at    which  99%   of  the   free-stream  den- 
sity  is   recovered.      Finally,    the  sheath    can   be   self-generat- 
ing  from   a    consistent   set   of  equations   without   an   overall 
energy    equation. 

The    resulting   electric    fields   at   the   free-stream   boun- 
dary  do    not   quite   attain  the  desired    2.7x10s   7/m    presented 
in   the    problem   formulation.      These   results   fall   short   by  the 
same   factors  required    for   the   adjustment    of   the  electron 
temperature   field   to   meet   the   stated   free-stream   boundary 
requirements.      It    must   be   pointed   out    that   the   electric 
field  is  the  required   boundary   condition.      However,    this 
boundary  condition   could  be    made  to   yield    no  reasonable 
results,   and   so   the   chosen    boundary   condition  is   ^  =    0. 

The   resulting   electric    fields  at   the    active   site   in  each 
of   Cases   I,    II,    or   III   was    of  two   orders    of   magnitude   higher 
than   that   of  the   free-stream.      These   electric  fields    are 
below   that    needed   for    breakdown   of   Nitrogen   at   one   atmo- 
sphere,   which  is    determined   as,    [Ref.    18]; 

2fr    =    (E/p)  p   =    (350   V/cm-torr)  x(760   torr)    =   2.7x107   V/m 
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The  requirement  of  strictly  matched  currents  produces 
steep  derivatives  and  inflections  of  the  density  profiles  in 
the  vicinity  of  the  anode.   A  noticeable  "spixe"  is  evident 
in  Pig.  12- j  for  example.   Case  II  and  III  anode  electron 
densities  were  established  on  more  aesthetic  grounds,  where 
smoother  appearing  density  profiles  were  generated.   In 
these  cases  the  anode  current  exceeded  the  free  stream  cur- 
rent by  only  a  factor  of  3. 

The  nonactive  regions  of  the  electrode  (wall)  far  from 
the  active  site  reveal  a  small  positive  space  charge  den- 
sity amounting  to  approximately  2%   of  the  free-stream  den- 
sity.  The  negative  space  charged  sheath  is  46%  to  56%  of 
the  free-stream  density. 

Previous  work  [ Ref .  12]  of  a  similar  problem  reported  no 
reasonable  results  for  the  case  of  a  catalytic  wall 
(n^=n-=0).   The  results  presented  in  this  study  demonstrate 
wall  conditions  of  nc- =  0.  and  n^  <  1 .   Here,  a  noticeable 
positive  space  charge  coexists  with  the  large  negative  space 
charged  sheath,  without  large  inflections  in  either  the  vol- 
tage potential  or  densities. 

The  Case  III  study  has  both  electron  and  ion  wall  condi- 
tions floated.   The  shape  of  these  density  distributions  is 
quite  symmetrical  about  the  anode  site  in  all  directions. 
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One  needs  little  imagination  to  visualize  that  the  density 
profiles  can  be  "wrapped"  around  the  anode  completely,  in 
effect  treating  the  anode  as  a  whisker  or  protrusion  into 
the  plasma.  However,  the  imagination  must  be  tempered  by 
the  fact  that  it  can  only  be  a  two-dimensional  protrusion 
(e.g.,    a   "blade") . 

The   effect   of   the   net    production   term   has  been   previ- 
ously  discussed   as   one  of   the   most   important  results   of  this 
study.      The  term   is    necessary    for   inclusion  in   the    problem 
in   order    (1)    to   attain  wall   electron   density  values   of   less 
than   the   free-stream,    (2)    to  produce   the   boundary-layer 
behavior  of   the   species   densities, and    (3)    to  obtain   con- 
verged  solutions   with   voltages  greater   than    10   volts   or   com- 
putational  node   spacings   of   any   order   less   than   that   of  the 
sheath   length. 

C.       NUMERICAL   CONCLUSIONS 

The   Modified   Newton-Raphson    (MNR)     method,    presented  to 
solve   the   particular   set   of   coupled   non-linear   partial   dif- 
ferential equations   in  this    study,    is   a   general   method   that 
can   be   applied  to    more  sophisticated   systems.      The   defi- 
ciency  of   the  quasi-Jacobian  method    (discussed   in   Sec.    III.C 
and  App.    C.)    in   dealing  with  non-linear   terms  is   circum- 
vented  by   an   alternative  linearizing   scheme.      Convergence  is 
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rapid;  however,  the  starting  solution  must  not  exceed  10 
volts,  nust  have  zero  net  production,  and  must  have  a  cons- 
tant electron  temperature  during  the  first  2-3  iterations. 
After  the  solutions  take  their  initial  form  at  ten  volts, 
the  potential  field  is  "amplified"  and  reinstated  as  a  new 
starting  solution.   The  final  stages  of  the  converging  solu- 
tion are  carried  out  under  the  full  influences  of  (1)  float- 
ing wall  densities,  (2)  variable  electron  temperature,  (3) 
net  production  and  (4)  an  active  node  site  electron  density 
being  solved  for  through  current  matching. 

The  profiles  of  the  space  charge  and  net  production  dis- 
tributions demonstrate  erratic  behavior  at  the  first  few 
computational  nodes  from  the  active  site.   This  is  believed 
to  be  caused  by  the  linearization  of  the  governing  non-li- 
near equations  and  has  the  appearance  of  truncation  error 
when  observed  on  an  iteration  to  iteration  basis.   However, 
tne  remaining  regions  of  the  profiles  are  smooth. 

D.   RECOMMENDATIONS 

This  work  succeeded  in  solving  the  problem  in  a  two-di- 
mensional environment.   This  fact  necessarily  limits  the 
quantitative  results  since  the  current  density  decreases  in 
only  two  dimensions.   with  the  computational  power  of  recent 
computers,  it  should  prove  worthwhile  to  make  the  next 
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advance   as   a   three-dimensional    problem   rather  than   elaborate 
in   the   two-dimensional  case.      &    three-dimensional    scheme 
would   render   more   realistic   results  of   a  quantitative    nature 
because   active   spots   are  more    likely   to   exist  than   active 
lines  and   because   3-D   geometries  require   less  voltage    for  a 
given   current. 

Specific  recommendations  other   than  the   above    include; 

•  Account    for   the  ion   contribution   to   the   current 
boundary   condition    into   the   wall. 

•  Investigate  the  physical  implications  of   the    low 
electron   temperature  near   the   wall   far   from   the 
anode  site. 

•  Investigate  the  physical   implications  of   the   elec- 
tron  density   "spikes"   prevalent   near  the   anode  site. 

•  Make   the    wall   temperature   a   valid   boundary   condi- 
tion. 

•  Use   a   free-stream   boundary   condition  of 
E^>=2.7x105    V/m,    rather   than   zero   potential.      This 
would   introduce  desired  results    for    pumping    lasers. 

•  Include  a   more  sophisticated   electron   temperature 
model,    or   verify  the  present   one. 

•  Increase   the   computational   mesh   density. 

•  Study  the   general   application   limits  of   the    MNR 
method. 

•  Investigate   "n  "    in   the   role    of   making   one-dimen- 
sional  solutions  possible    (i.e.,    alter  the   reacting 
rates) . 

•  Investigate   whether    the   low    electron  temperatures 
and  associated  steep  derivatives   located   along  the 
wall   between   active   sites  have  any  physical    implica- 
tions,   or  result   merely   from   inadequate   problem   for- 
mulation. 

•  Consider    the   electron-beam   ionization   in   the   solu- 
tion to   produce  more  realistic  effects. 
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APPENDIX  A 
CHARACTERISTIC  PARAMETERS 
A.   SHEATH  LENGTH 

The  extent  of  the  sheath  is  one  of  the  most  important 
charcteristic  lengths  in  this  study.   It  is  within  the 
sheath  that  most  of  the  potential  drop  occurs  for  a  short 
discharge.   Fortunately,  the  sheath  length  can  be  estimated 
rather  easily  from  Poisson1 s  equation,  (Eg.  A1) 


Let; 


?V  =  f:(we-ru)  (i1) 


*••■*&■' -i£  '    *>—**       » 


where  "  \*m    is  the  characteristic  sheath  length;  the  sub- 
script "o"  indicates  the  value  at  the  electrode  and  "   "  the 
value  at  the  undisturbed  plasma.   After  gathering  the  dimen- 
sional quantities  on  one  side,  Eg.  A2  becomes; 

7*0  s   &S§eJ3.iWea"v>«\  (A3) 

The  characteristic  sheath  length  ■ fa,    defined  in  this 
fashion,  is  used  to  non-dimensionalize  the  "del"  operator  in 
all  directions.   That  is,  independent  of  the  dimensionality 
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of  the  problem,  the  sheath  is  effective  over  the  character- 
istic length  just  defined.   Now  an  estimate  of  the  size  of 
the  sheath  can  made  by  setting  the  grouping  of  dimensional 
quantities  in  Eg.  A3  to  be  of  the  order  one,  namely; 

where  "  h^jj    is  the  Debye  length  in  the  undisturbed  plasma. 

Anticipated  values  are  now  established  for  the  sheath 
length.   Take  for  example  a  35  volt  anode  potential,  and  a 
charge  species  density  of  1017  to  1019  m-3.   These  parame- 
ters define  a  sheath  length  of  1.4x10-*  to  1.4x10-s  meters. 

As  can  be  seen,  a  ten-fold  increase  in  voltage  would 
change  "Jcn  by  a  factor  Df  3,  rendering  the  nondimensionali- 
zation  of  the  problem  relatively  insensitive  to  the  chosen 
value  of  the  voltage  potential.   That  is  to  say  the  value  of 
the  sheath  characteristic  length  as  defined  may  be  consid- 
ered a  resonabie  estimate  for  the  charge  species  densities 
examplif ied. 

The  important  implication  here  is  that  extreme  care  must 
be  devoted  to  electrode  surface  preparation  to  preclude 
whiskers  or  flaws  of  a  size  comparable  to  " ^$".   These  pro- 
trusions would  produce  a  focusing  of  the  electric  potential 
to  create  extremely  high  electric  fields,  which  become  sites 
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for  initiation  of  breakdown.  Protrusion  growth  has  been 
observed  to  form  in  a  time  on  the  order  of  milliseconds. 
Many  conditions  are  favorable  to  protrusion  growth,  depen- 
dent on  the  pre-breakdown  and  breakdown  current  profile, 
electrode  composition  (especially  in  the  case  of  a  ther- 
mionic cathode) ,  surface  purity  and  surface  preparation 
[Hef.  21].   Some  possibilities  axe; 

•  Electrode  plastic  flow  or  meltinq  in  an  electric 
field  due  to  micro-particle  bombardment. 

•  Sticking  of  microparticles. 

•  Electron  beam  induced  protrusion  growth. 

i 

B.   DIFFUSION  LENGTH 

If  diffusion  is  a  significant  factor  in  this  study  then 
its  contribution  to  the  current  would  be; 


Jv  =  D*  e 


an*  -  0,  e^i 


■*-??- 1<5 


or,    since   D-    >  D«,    and  using  Einstein1 s   relation,    or 
D/j4€=kTc/e; 

making    the    equation    unit   free    (through  fractional   analysis) ; 

JL-       .    £rvr    _  ottDefrM  (A5) 


e 


Vt  **+         9  x 
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If   the   right  side   of   Eq.    A5    is    of   the    order  one,    (the  nondi- 
mensional   density   gradient    is  a    maximum   of    1)    then   "LM    may 
be   estimated   as; 

L    a  eDen-    /J  «.  3.7x10-*   m 
for   the    values; 

ne   =    10la   m~3 
De   =   0.0665    m2/(V«sec) 

Jc  =^»nft»^^  s</lf«e«ae«n0«(E/n)  «s    2.88x103    A-m~3 
Since   ^>  Lr    then  any   role   attributed   to  diffusion 

should    be   visible   within  the  sheath.      The    E/n   within   the 

thin   sheath   is   of   course  greater  than   the    E/n  at   the   free 

stream    (10*0    v«m2)#    which    decreases   the   local  value    for    "L" 

stiil  further. 

Another   approach   is  to   assume   a   one-dimensional  solution 

to   the   charged  species   conservation   equations,    which    can 

only   be    done  outside   of  the   sheath,    i.e.,    in   the   ambipolar 

region.      The   ion    conservation  equation  may    be  written; 

n   =  7#J  =   e7«C7nW  -   Dv'n]  (A 6) 

Since   it   can   be   shown    (next    section)     that   the  production 

term  can   be   approximated  as; 

n  aiP^n   -OCn2 

then  for  the  stationary  condition  of  no  net  production; 

?);  -  ot  N 

Hondimensionalizing   Sq.    A6,    it    can   be   shown   that   the 
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characteristic  diffusion  length  ML"  is; 

L  =  J  b/tfN  =  2.x10-«  meters 
which  is  of  the  same  order  of  magnitude  as  the  characteris- 
tic sheath  length. 

C.   IONIZATION/RECOMBINATION  LENGTH 

Analysis  of  the  production  term  (Eg.  4)  reveals  that  for 
Nitrogen  the  three- body  (electron)  collisional  recombination 
rate  (tf-J  is  small  in  relation  to  to  the  other  terms,  (see 
App  B.) .   In  any  case  the  net  production  term  is; 

nc  =  *}n€-<*;nent.  -  ajn^n,  n0  (A7) 

In  the  steady  state  the  net  production  term  is  zero,  and 
as  is  done  in  Appendix  B,  the  relation  becomes; 

ne  =    (#,♦  0^n#)  •  (&-   n^)  «n^ 

M 

substituting   the    value  of    0.17x10-* l    a3/sec,    and    (Er  -    ru )  a 

1018,    we  can  estimate   the   characteristic   production    time  to 
be    Tp  =    5.88x10-7    sec. 

To   establish    a   meaningful  characteristic   length,    the 
average   thermal   velocity   will   be   used   since   it   represents 
the   principal  transport   velocity   of   the   medium. 

Jp  =  Tf>*vl    =    5.88x10-*«3.2x10*    =    1.88x10-*    m 
The   characteristic   length  for   net    production   is    a   factor 
of    four    larger   than   the  sheath    length    (A<=   4.4x10~s    m) .      The 
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conclusion  which  can  be  made  from  this  (fractional)  analysis 
is  that  the  net  production  is  significant  within  the  overall 
region  of  interest,  including  the  sheath  itself. 

D.   FROZEN  PLOW 

When   a    process   relaxation  time   "tv"    is   long   compared 
with   the  time  that  a    fluid    particle   spends  in  a   given 
region,    the   flow   is    said  to    be   "frozen"   in   this   region.      By 
this  term  one   means   that  a    reaction  process   may   be 
neglected,    as   far    as    influencing  the    flow   is  concerned,    and 
the   species   involved   treated  as   chemically   inert,    [Ref.    22]. 
In  the   opposite   limit,    where  "Tp"  is    very   small   relative  to 
nTrn   the  characteristic   flow   time,    the   flow  is   said   to   be   in 
"equilibrium."      For   the   recombination   rates   given    in   this 
section,    "Trn   =   6.x10~r   and    nt-,n   -   4.0x10~7   seconds.      Here 
flow   time   is   defined   as   the   sheath   characteristic   length 
divided    by   the   drift    velocty.      Thus   the   flow   is   not   "frozen" 
as    far    as   chemical  reaction,   but   we    must   consider   that   the 
electron   and  ion   species   are  not   in   thermal   equlibrium. 
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APPENDIX  B 
DEVELOPMENT  OF  THE  PRODUCTION  TERM 
The  net  production  term  is  defined  as  the  difference 
between  ionization  and  recombination  rates.   Whenever  a  neu- 
tral atom  or  molecule  loses  one  or  more  of  its  electrons, 
ionization  has  ocurred.   Recombination  is  the  process  wher- 
ein positive  ions  and  electrons  combine  to  form  neutral 
atoms  or  molecules.   It  is  simply  reverse  ionization.   If 
the  net  production  is  zero  then  ionization  and  recombination 
are  equal.   The  net  production  term  is  the  right  side  of 
2qs.  2  and  3,  the  expression  for  which  is  developed  here. 
This  work  considers  only  three  types  of  recombination 
(1) Dissociative  -  the  electron  recombines  with  a  molecu- 
lar ion  (AB)  +  ,  and  the  recombination  energy  goes 
into  dissociating  the  molecule  and  increasing  the 
kinetic  energy  of  the  resulting  products.   The  pro- 
cess is  described  as; 

e  ♦  (A3)  +  — *r  A  ♦  B 
(2)  Three- body  (electron)  -  the  energy  released  in  the 
recombination  process  is  carried  off  by  a  third  body 
which  is  an  electron.   The  process  is  described  as; 
e  ♦  A*  *  e  — 9>   A  *  e 
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(3)     Three-body    (heavy^    -  the  recombination   energy   is 

carried   off   by  a   third   body   which    is  a   heavy    parti- 
cle  such    as    an   atom    or    molecule.      The  process   is 
described  as; 

e   ♦   k*    +   B   — ¥    A    ♦   B 
This   study   does   not  consider  other  recombination    mechan- 
isms  such   as  radiative  or   dielectric   recombination,    since 
they  are   generally   of  low   significance  in    relation    to   the 
just   described   processes.      The    ions  are   limited  to   singly 
positive  ions.      No  account    is  made   for  intermediate   or 
excited   states   as   this  requires   a   detailed   balancing   which 
is   beyond  the  scope   of  this    work.      Additionally,    no   account 
is   made   of    attachment. 

Rewriting   Eg.    4    for   simplicity   of   discussion; 

&e   =    ^ne  '^i^i    "  *2ne2nc    ~ ** ne n<:  no  <B1> 

where   each    "^•',    refers  to   the  recombination   coefficients  for 

dissociation,    three-body    (electron)    and   three-body    (heavy)  , 
respectively. 

The    recombination   coefficients   depend   upon   the   ionic 
species    involved.      Some  general  but   representative    values 
[Ref .    22  ]  follow; 

tf,  =    10~13    mVsec 

o^x=   9x1  ()-*♦    m3/sec 

£.   =   2x10-13    mVsec    (air) 

92 


Hinnov   and   Hirschberg   [  Ref .    23]   developed  a   relation  for 
three- body     (electron)    recombination   coefficients   as   a   func- 
tion  of    electron    temperature    (Tfc )    and    density ,    valid    for 
energies  less  than  0.25   e7. 

Ux~    1.09x10-*on    t^-V2    mVsec 
where   the  electron   temperature    is   °K.      Gurevich   and   Pitaen- 
ski  [Ref.    24]  extended  this    development   to   obtain; 

*i=    2.  13x1 0-2° neTe-VsinwV   m3/sec 
where; 

lnwi=    1.24x10^     (Te3/ntf)-s 
At    1.0   eV  these  coefficients   are  5. 6x10-2'    and    1.05x10-*° 
i3/sec   respectively.      A  simplification   of   the  later   expres- 
sion  for    #\  is   valid   from   energies   of    0.1    to    10.0   eV   and 
n&=10-i«   m-3; 

<*z=0.0  32xT-'»-  3 
Brown  [Ref.    18:    p. 197]    gives  a    value    for  the   dissocia- 
tive  recombination   coefficient; 

<*,    =     1.  5X10-12     m3/S 

The  ionization  coefficient  can  be  determined  from  con- 
sidering that  there  is  no   net  production  in  the  free  stream, 
i.  e.  ,    far  from  the  electrode.   Here; 

0  =*An  -<*,N2  -a^N3  -*jN2n#  (B2) 
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where  "NH  is  the  stationary  state.   Solving  for    and  sub- 
stituting into  Eg  31  gives; 

ne=  n£(N-n-  )  (*,*0^no)  *oCln9,T$*/**  (Nz-n^nj  )       (B3) 
where; 

*,  =    1.5x10-12  (from   Ref.    18) 

4^=    1.09x10-20  (from    Hef.    23) 

0^no=   2.0x10-13  (from    Ref.    22) 

The    net   production   term  is    the   relation,    suitably   non- 
diaensionalized,    that   appears  in   the   computer  subprogram 
"RTSIDE. " 
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APPENDIX  C 
NUMERICAL  SOLUTION  METHOD 

A.   GENERAL  DISCUSSION 

There  are  two  basic  approaches  used  in  solving  a  nonli- 
near system  of  equations,  such  as  possed  in  this  study.   A 
given  method  can  usually  bs  classified  as  one  of  the  ascent 
(or  descent)  or  more  commonly  refered  to  as  gradient  meth- 
ods, or  as  one  of  the  Newton  methods.   The  method  applied  in 
this  study  and  by  Dolson  [Ref.  12]  is  one  of  the  Newton 
methods.   A  brief  review  of  both  methods  will  be  presented 
to  facilitate  development  of  the  particular  method  employed 
in  this  study. 

Let  the  nonlinear  system  of  equations  be  written  as; 

f/  (  ^  |  r  X  j  g  •  •  •  #  Xj)  )   =  0 

.  .  .  (C1) 

f n  ( x  f  *  x2  »  •  •  •  '  xf) )  = 

or  in  vector  form  for  convenience; 

f  (?)  =  0  (C2) 

In  the  descent  method  the  problem  is  possed  in  such  a 
way  so  as  to  form  a  function  which  is  zero  at  any  solution 
and  positive  otherwise.   A  typical  such  function  would  be; 
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n 
0(xi,x2,...xn)  =Xl|fi|2  <c3> 

so  that  the  problem  is  reduced  to  finding  a  minimum  of  0, 
which  also  happens  to  be  the  value  zero.  In  solving  prob- 
lems by  this  method,  an  initial  guess  to  the  solution  must 
be  provided,  call  it  x<»>.  Then  a  direction  d  is  determined 
so  as  to  point  to  a  better  solution,  x<2>.  This  direction 
is  often  chosen  to  be  the  gradient  of  0(x).  The  operation 
-grad  0(x)  is  the  direction  in  which  0(x)  decreases  most, 
and  is  locally  the  best  direction  in  which  to  reduce   0(x) . 

d  =  -grad  0(x) 
so  that ; 

x<2>  =  x<i)  -  a«grad  0(x<*>) 
where  "a"  is  chosen  to  be  an  appropriate  step  size.   A  prob- 
len  is  encountered  at  local  minima  of  0(x) ,  which  needs  to 
be  avoided.   One  must  also  use  the  step  "a"  in  such  a  manner 
as  to  control  the  acceleration  along  slowly  or  rapidly  vary- 
ing 0(x).   It  may  be  that  by  the  nature  of  0(x)  ,  the  gra- 
dient may  have  to  be  approximated  by  finite  differences,  an 
approach  which  offers  a  flexible  recourse  even  when  not 
reguired. 


96 


The    Newton    method    involves    approximating   the    nonlinear 
problem    by    a   suitable   linear   problem    involving   derivatives 
of    the    nonlinear    equations    evaluated   at    x<l>.      The   linear- 
ized   problem   is   simply  a   set   of    linear   algebraic   equations/ 
the   solution   to   which    furnishes  an   increment   that    when   added 
to    the    solution    x<l>    will    yield    x<2>.      This   method    is   a    gen- 
eralization  to   n-dimensions    of    Newton*s   method   for    evaluat- 
ing  the    zeros  of   a   function. 

Each    function    f±  (x)    is    approximated   by    the    initial   two 
terms   of   a   Taylor   series   expansion   about    the   solution    vector 
x<  i>  ; 

fi(x)    =    fi(x<i>)    ♦      J(x<i>>«[x-    x<*>]  (C4) 

where   "J  "    is   the    familiar    Jacobian   matrix, 


4-4  a*j 


*")  <C5) 


j  =  l 
The   right-hand-side   of   equation   C<*    is   the    linear   vector 

function    (tangent  hyperplane)    of   x   which    best   approximates 
the  nonlinear    function   f   evaluated   at    x< l > .      He   set    this 
right-hand-side    equal   to   zero   and   solve    for    x,    which    is  an 
improvement    from    x< l > .      He    call    this    improved   solution    vec- 
tor   x< 2 >  .    and    from    equation    C4    we    j*t; 

x<2)    =    K(i»    -       J(x<  l>)  "l»f  (x)  (C6) 
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The  Jacobian  need  not  be  known  very  accurately,  and  often  is 
left  unchanged  for  several  iterations  since  its  evaluation 
often  entails  penalizing  computer  time.   Thus  the  method 
(slightly  modified  as  described  briefly  in  section  III.B) 
used  by  Dolson  [Ref.  12 ],  is  referred  to  in  many  texts  as 
the  Newton  method  or   as  the  Newton-Raphson  method  (Kefs.  25, 
26]. 

The  Jacobian  matrix  of  coefficients  for  elliptic   prob- 
lems is  typically  a  sparse  banded  matrix,  to  which  there  are 
many  computer  subroutines  available  for  solution. 

B.   MODIFIED  NEWTON-RAPHSON  METHOD 

The  particular  solution  technique  employed  in  this  study 
merits  a  detailed  discussion,  since  it  is  not  common  in  the 
literature.   Bailey  and  Touryan  (Ref.  27]  use  a  similar 
method  and  provide  some  description  of  its  development.   The 
method  is  similar  to  other  Newton  methods  however,  the  Jaco- 
bian matrix  is  replaced  by  a  matrix  of  coefficients  devel- 
oped from  a  sequence  of  linear  differential  equations  which 
approximate  the  system  of  nonlinear  equations.   This  method 
will  be  referred  to  as  a  modified  Newton-Raphson  method  or 
MNR  method  for  short,  and  will  be  presented  in  general  using 
a  specific  example. 
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Take  Poisson1 s  equation  given  earlier  and  repeated  here 
for  simplicity; 

The  known  solution  variables  are; 

£<*>,  ne<t  >,  and  nc  < * > 
In  general  the  first  stej>  is  linearization  of  the  equa- 
tion.  In  this  case  Eg.  21    is  linear,  but  it  will  serve  as 
an  example.   ie  assume  that  the  suceeding  solution  variables 
at  iteration-  (2) ,    namely  z<2>,  can  be  approximated  by  the 
solution  z<*>  plus  some  small  change  ^z<l>. 

4Z<2>  =  Z<1)  +     ZC1)  (C8) 

where; 

It  remains  to  determine  4z(1).   Substitution  of  this 
relationship  into  the  Poisson  Equation  is  the  second  step; 

vl<»a) ,  vVV  vW>«  c,K'-  *rt)  *  c,  (W- W) 

rearranged; 

7V'-C,(A<' -  A«iu>)  =  -  [vV°-  c,  (n«°-Oj   (C9) 
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The  right  hand  side  of  Eg.  C9  is  a  function  of  the  known 
solution  vector  z<»>.   We  are  applying  this  relationship  at 
each  computational  node.   A  second  order  partial  differen- 
tial equation  can  be  approximated  by  a  second  order  differ- 
ence equation  or  "3-point  formula,"  so  Eg.  C9  will  become  a 
system  of  linear  algebraic  relationships  involving  near 
neighboring  values  of  the  solution  vector  z{1>.   The  left 
side  of  Eg.  C9  is  composed  of  a  system  of  coefficients 
expressible  in  terms  of  z<i>  times  the  unknown  (perturba- 
tion) vector  £z<*>.   So  we  now  have; 

A(Z<  l>)*Az<  i>  =  -F(Z<*>)  (C10) 

Here  the    matrix   of  coefficients   "A"    is   a  sparse   matrix 
with   a   relatively    narrow   band   of  non-zero   elements   about   the 
main   diagonal,    provided  that   the  "natural"    order  of   the 
resulting   linearized   equations   is  preserved.      There  are 
standard  computer   routines    available   which   take   advantage  of 
this   banded   structure,    such   as   LEQT1B,    an   IHSL   routine. 

It   follows   that  since    4z<*>   =   z<z>    -    z<*>; 

z<2)    =  z<i)   -  r  A(Zct>)  ]-i«F(z<")  (C11) 
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For  the  case  of  a  nonlinear  term  such  as; 

there  is  no  additional  difficulty.   when  linearized  this 
term  becomes; 


fte.^  tAJIe^^  Y)*d£*  *    ^*e<^  (C12) 

dX      <^x      c5x        dx 


He  neglect  high  order  terms  created  by  the  inner  products  of 
the  elements  of  Az<*>,  such  as  the  last  term  in  relation 
C12.   The  first  term  would  be  a  part  of  the  function  F(z<l>) 
in  the  right  side  of  Equation  C10  ,  and  the  remaining  two 
terms  would  be  included  as  part  of  the  expression: 

A  (z<i  >)»Az<n. 
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APPENDIX  D 
DESCRIPTION  OP  COMPUTER  PROGRAMS 

This  appendix  serves  as  a  user  oriented  guide  to  the 
computer  programs  listed  in  App.  E.   The  programs  are  listed 
in  the  order  used  to  generate  solutions  to  the  problem  for- 
mulated in  this  report.   The  output  consists  of  printed 
results  for  each  iteration  and  punched  output  to  provide  a 
permanent  file  copy,  which  is  to  be  used  as  input  to  the 
three  plotting  programs  also  listed  in  App.  E.   These  plot- 
ting programs  use  the  Vecsatec  plotter  and  associated  soft- 
ware, and  represent  the  only  facility  not  commonly  available 
elsewhere.   The  main  program  SHTH  as  presented,  is  config- 
ured to  also  produce  output  at  a  monitor  terminal  as  used  on 
CP/CMS  (I3M  360)  or  MVS  (IBM  370).   The  output  in  summary  is 
via  three  devices; 

DEVICE  4  -  typically  a  CRT 

DEVICE  6  -  typically  a  printer 

DEVICE  7  -  typically  a  punch  or  disk  output 

These  can  be  manipulated  in  convenient  configurations  to 

give  a  full  or  partial  (sampled)  display  of  the  solution  at 

each  iteration,  or  merely  to  display  convergence  information 

in  a  simple  format.  - 
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The  SHTH  program  is  the  main  program  which  together  with 

its  9  subroutines  and  1  IMSL  subroutine  solves  the  non-di- 

aensionalized  forms  of  Eqs.  1  through  4.   The  input  to  the 

program  consists  of  control  cards  and  data  cards,  namely; 

1  control  card 

3  or  104  initial  solution  data  cards 

1  control  card 

The  control  card  is  described  here  in  terms  of  the  Fortran 

variable  names  thereon. 

ITHIH  -  is  the  iteration  counter  (ITER),  starting  value. 
If  ITHIH  ■  1,  then  the  program  will  read  in  the 
next  3  cards  as  data  which  contain  the  short  form 
starter  solutions.   If  ITHIH  is  greater  than  1, 
then  the  program  will  read  in  the  next  104  cards  as 
data,  representing  a  detailed  starting  solution. 
This  detailed  starting  solution  is  typically  the 
output  of  a  previous  run  of  the  program.   If  ITHIH 
=  0.0,  then  the  program  will  not  read  in  any  data* 
being  thus  informed  that  the  problem  is  over. 

ITHAX  -  is  the  last  iteration  value.  The  Droaram  ceases 
iterations  at  this  point  and  outputs* the  solutions 
and  auxiliary  information. 

ITHDT  -  enables  the  employment  of  the  net  production  term 
(BTSIDE)  when  ITER  is  greater  than  ITNDT,  otherwise 
the  net  production  is  set  to  zero. 

ITEHP  -  for  ITER  less  than  ITEMP,  the  electron  tempera- 
ture is  isothermal,  otherwise  it  is  a  function  of 
E/n.   (see  SHTH  statements  899-990). 

ITBC  -  for  ITER  greater  than  5  and  less  than  ITBC,  the 
wall  boundary  conditions  for  the  species  is 
floated. 

LADJ  -  1   electron  temperature  is  calculated 
-  0   electron  temperature  is  isothermal. 

LPOI  -  1   output  is  punched  on  device  7.   -  0   output  is 
not  punched. 

H  -  sets  the  mesh  spacing  equal  to  H*LAHS  (sheath 
length) . 

FKTR  -  is  the  adjustment  factor  to  the  temperature  equa- 
tion (see  Sec.  IV)  and  is  used  when  LADJ  =1  and 
ITER  is  less  than  ITEHP.   The  adjustment  factor  is 
not  used  when  LADJ  is  set  to  zero.   When  the  temp- 
erature is  coupled  and  LADJ  =  1 ,  then  the  program 
calculates  its  own  adjustment  factor. 

NAHE  -  is  reserved  for  user  labelling  and  identification 
which  is  used  at  output. 
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The  solution  set  is; 

PH  voltage  potential  (volts)  . 

HE  electron  density  (10*8  m*3) 

HI   ion  density  (1018  m-3) 

TH  nondimensional  electron  temperature 
The  short  form  starter  solution  is  illustrated  in  App. 
2.   It  consists  of  PH,  HE  and  HI  values  at  the  anode,  wall 
and  equilibrium  boundaries.   The  long  form  starting  solution 
consists  of  PH,  HE,  HI,  and  TH  distributions,  each  one 
represented  by  25  data  cards  and  a  label  card.   Each  card 
has  7  values,  (see  fomat-803) . 

A  typical  iteration  secies  would  start  out  for  about  ten 
iterations  from  the  short  form  data  with  ten  volts  at  the 
anode  and  other  appropriate  boundary  conditions.   The  ini- 
tial solution  converges  very  well  for  H  =  0.5,  the  net  pro- 
duction term  coupled  after  2  or  3  iterations,  and  the  wall 
conditions  not  floated.   This  output  is  then  used  as  input. 
The  voltage  is  amplified  to  say  25  or  35  volts  and  a  new 
control  card  is  structured  to  couple  the  temperature  or  wha- 
tever else  is  desired.   The  H  value  may  also  be  changed  to 
reflect  a  smaller  mesh  spacing.   It  must  be  kept  in  mind 
that  too  many  changes  in  one  solution  run  can  be  catas- 
trophic. 
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Each  iteration  of  the  solution  by  the  SHTH  program 
accomplishes  the  following; 

(1)  floats  the  NE,  NI  wall  values 

(2)  sweeps  through  the  A-matrix  and  C-vector  setting 
up  the  equations  that  satisfy  the  boundary  condi- 
tions and  the  internal  mesh  conditions 

(3)  couples  the  electron  temperature  equation 

(4)  couples  the  net  production  term 

(5)  calculates  the  anode  electron  density  required  to 
satisfy  current  matching  btween  the  anode  and  the 
free-stream. 

The  system  of  equations  has  the  form  presented  in  Eq. 
C10  in  App.  C,  and  is  solved  by  use  of  the  iflSL  routine 
LEQT1B.   The  solution  vector  is  then  updated  from  the  return 
information  contained  in  vector  C. 

There  are  two  solution  convergence  monitors.   The  first 
is  the  sum  of  the  squares  of  the  differences  (spectral  norm) 
of  the  -F(z<x))  values  from  Eg.  C10.   The  other  is  the  spec- 
tral norm  of  the  Azci)  upon  return  from  LEQT13.   Both  of 
these  are  available  at  the  video  terminal  (device  4)  and 
printer  (device  6,) . 

After  solution  update,  the  elctron  temperature  is  calcu- 
lated and  adjusted  per  the  control  card  options.   It  should 
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be  pointed  out  that  on  the  last  iteration,  the  temperature 
distribution  is  automatically  calculated  for  output,  but  is 
not  used  in  the  scheme,  unless  otherwise  coupled  per  the 
control  card. 

The  matching  of  the  anode  current  to  the  free-stream 
current  is  accomplished  by  solving  for  the  anode  electron 
density  using  subroutine  ANODE.   The  call  to  this  subprogram 
will  require  a  change  to  the  SHTH  program  as  it  is  presently 
overridden  with  a  comment  card,  (prior  to  statement  9000)  . 

The  final  outputs  are  straight-forward,   subroutines 
Jf ALL  and  JO0LE  are  called  in  this  phase  of  the  program,  in 
order  to  compute  the  current  into  the  wall  and  the  Joule 
heat  distribution.   The  electric  field,  net  production  and 
space  charge  density  distributions  are  also  calculated  and 
provided  as  punched  and  printed  output. 

The  punched  output  is  used  as  input  for  the  STHEAH,  PICT 
and  PLOT  programs.   These  programs  use  the  Versatec  plotter 
to  provide  a  presentation  of  the  current  stream  lines,  an 
oblique  perspectives  of  the  solutions  and  plots  of  the  solu- 
tions as  they  vary  normal  to  the  wall  extending  from  the 
anode  site.   The  details  of  these  programs  are  combersome, 
and  it  is  wise  to  leave  them  to  the  user  for  inspection. 
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The  first  two  distributions  in  the  punched  output  are 
used  in  the  STB2AM  program,  requiring  that  the  nonpertinent 
data  be  removed.   The  next  8  distributions  are  used  for  the 
PLOT  and  PICT  programs.   The  PICT  program  has  an  option  card 
filled  with  ones  or  zeros  in  the  first  eight  positions. 
This  option  card  enables  the  user  to  elect  which  distribu- 
tions to  plot.   The  plots  provided  in  this  report  as  Figs. 
12  through  14,  were  produced  using  the  PICT  and  PLOT  pro- 
grams.  The  programs  are  simple  and  speak  for  themselves. 


107 


APPENDIX  E 
COMPUTER  PROGRAM  LISTINGS 
This  appendix  contains  a  listing  of  the  programs  used  in 
the  solution  of  the  equations  presented  in  this  work 
(namely:  Egs.  9,10  and  11),  and  the  auxiliary  programs  that 
provided  the  plots  for  the  report. 


SHTH  (MAIN  PROGRAM) 109 

DELOP 1  27 

GRDNT 129 

TEMP 131 

RTSIDE 132 

PATCH 1  33 

ANODE 134 

J  WALL 136 

JOOLE 137 

OUTPUT " 1  38 

PLOT 1  39 

PICT 141 

STREAM 144 

INITIAL  INPUT  (sample) 150 
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